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I. INTRODUCTION 


A wide variety of methods are available for determining the pressure 
field from a point source in a stratified ocean of constant depth. 
Figure 1.1 Shows that, by using various transforms, it iS possible to 
derive expressions describing the pressure field in terms of the linear 
wave equation, the Helmholtz equation and, finally, the field integral. 
Methods associated with these expresssions vary, and each has relative 
advantages and disadvantages in terms of speed, accuracy, and ability to 
deal with complex situations. Unfortunately, no single method is supe- 
rior in all respects and there must be trade-offs. In some situations, 
accuracy might be sacrificed for speed, or the ability to deal with 
unusual boundaries may override the importance of speed or accuracy. 

Figure 1.2 shows the common methods of solving the field integral to 
determine the pressure. The Normal Mode Solution, the Method of Stee- 
pest Descent, The Multiple Scattering Method and Direct Numerical Integ- 
ration are "exact" methods. The Fast Field Program (FFP), a method of 
numerical integration, must be considered an approximation, not only 
because it approximates a Hankel function by an exponential, but alse 
because of its reliance on a finite number of samples. Since Stationary 
Phase normally requires approximations, then so does its derivative, the 
Method of Images. The Wentzel-Kramers-Brillouin (WKB) Integral normally 
involves approximations so both Ray Theory and the WKB Mode Equation 


will seldom be exact. It is possible to find exact solutions using the 


Normal Mode Technique, but in this paper, as will be seen in Chapter II, 
practical assumptions relegate the solutions to the approximate cate- 
gory. This thesis is limited to consideration of the Normal Mode Solu- 
tion, Which is computed by the program EXACT, and a Fast Field Program, 
computer program FFP. 

Chapter II derives the field integral, then develops the Normal Mode 
Solution, and describes the Fast Field Program, FFP. The basic differ- 
ence between the methods involves the requirement of the Normal Mode 
Solution for an accurate estimation of the horizontal wave numbers for 
each of the modes. The FFP does not calculate the modal wave numbers 
DMiieecomputes the pressure contributions for a large number of pre- 
defined horizontal wave numbers. An expression to describe losses and 
phase shift in terms of the impedance and sound speed ratios between the 
bottom and the water is also derived in Chapter II. 

Chapter Ill describes the programs EXACT and FFP. EXACT solves ee 
field integral by first finding good approximations of the mode wave 
numbers using a matrix eigenvalue ete and then refining them. FFP 
Calculates the pressure field at a discrete number of horizontal wave 
numbers without regard to modal values. Both programs use the vector 
eum OL the components of pressure from the constituent wave numbers to 
calculate the propagation loss between the source and the receiver, 
taking into account absorption in the water, losses to the bottom, and 
phase changes at the surface and bottom. 

Chapter IV describes the testing procedures and the results. Within 
this chapter the testing criteria are defined and the adaptive nature of 


the development process is discussed. Test results provided some very 


instructive lessons which gave an insight into how normal modes can 
describe pressure variations in both the vertical and horizontal. 
Chapter V discusses the results and some of the weaknesses of the 


program. 


LO 


Exact wave Approx linate Geometrical 


SOUuULLoOns Solutions Acoustics 
Numerical 
a Solution 


Beste Phyics 
Cons of Mass 





—_————————— 
Method of Ray 
ee ee ee ee 
Characteristics Tneor 





Linear 





Lin Wave Eqn 








} Numerica 

>| Equation >| Theory 
[ Numerical 
Finite 
Blement 
Model 











Parabolic 
Equation 


Fourier Transform 
(t > w) 


| 


Helmnoltz Eqn 
V2P = K*P = F(w) 






Separation of 
Variables 
Pir ,z)=R(r)Z(z) 







! WKB | Ray 
> Approximation —» Theory 


| PETEMAG 302] 
ons 


SeeOOLiti 











ee 
Normal tiode 
Sei Ueren 


Couplea Modes | | Adiabatic | 
Pato (rju (r,2) IL»! Mode | 


So. UL oOns.. 








Hankel Transform 
(r+ &) 





Normal Mode 
Solution 





Field Integral 
Forme Source 


Eeeruy (1) ; 
a Ho (arjede 









Direct Numerical 
Solu oner fr 





PeetGe wet Methods of Solving £for the Pressure Field. 


2 





Field Integral 
(ir iy eevee 


W 





Direct 
Numerical 


Integration 





























steepest waveguide Normal | 
Descent Solution Mode 
d _ yee | CMULtLI ple Solutiam 
qe exponent aa Scattering) = 0 
2 1 | | | 










WKB Phase 
Integral 


Stationary 
Phase 


















1 
oe © 





Ordinary 
Ray Theory 






Figure 1.2. Solutions of the Wave Equation 
Derived from the Field Integral. 


12 


i GOR 


In this discussion of normal modes, frequency will be very low so 
the absorptive losses due to the water will be negligible, except at 
long ranges. However, absorption in the bottom will often be very high, 
especially for the higher modes. The theory was used to derive formulae 
for the lossless case; then absorption was introduced as an imaginary 
component of the horiztonal wave number. Although absorption in the 
water is low for low frequencies, it proved expedient to input abnor- 
mally high abdsorption rates into the test programs. Ey doing this. it 
was possible to verify that the complex numbers were producing consis~ 
ene values. 

The theoretical approach involves a series of transforms, some of 
which require approximations for simplification. These transforms 
Bese tvely Cake the wave equation from the time domain to the fre- 
quency domain, from three dimensions to two, from range dependent to 
Wave number dependent, and finally to a function that describes pressure 
as a function of the depths of the receiver and transmitter, the hori- 
zontal wave number and the horizontal range. Normal mode theory 
enables one to consider the total pressure as the sum of pressures from 
Many components. Each mode is stimulated by the transmitter, or sensed 
by the receiver, to a degree that depends on the positions of the 
transmitter and receiver in relation to the modal depth function curve. 


This curve is described in Chapter IV. 


ih 


In order to keep the program simple, several assumptions must be 
made; some limit the practicality of this particular srogram and comua 
be overcome in a more comprehensive routine, while others are of little 
consequence. 


The main assumptions, some of which will be discussed in detail 


later, are: 


1. Source is cw and monochromatic; it emits a continuous wave 
at a constant frequency; 

source is a point radiating uniformly in all directions; 
medium is homogeneous in all respects in the horizontal; 
bottom and surface are flat and parallel; 

surface is perfectly reflecting; 

bottom can be completely described by its impedance; 

all energy transmitted into the bottom is lost, and 

branch line integrals can be ignored. 


Ory AM Fw Ww 


A. THESP IED ite Chae 

An omnidirectional point source operating at frequency f(t) is posi- 
tioned at a depth z in a water mass that is uniformly homogeneous tin 
sound speed in. the horizontal plane and can be described in the vertical 
plane by the sound speed, c(z), at any depth. 

As long as the variations are linear, the pressure is glven by the 


acoustic wave equation: 


1 a 1 de i 
V“p Se ae f(tJéCs=s,) Zeal 


where: 


p = p(s,t) pressure as a function of position and time, 


f(t) = funetion describing the source amplitude, 
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e(z) = depth dependent sound speed, and 
S = position in Cartesian coordinates. 
A Fourier transform fGommtime to frequency, w, corresponds to con- 


Sidering the equation for distinct frequencies and yields: 


Veet 





Ww * Ete 
Bre P=" =F( a)'5(S—8,) eae 


where: 
P = P(S,w) a function of position and frequency, 
w = frequency in radians/sec, and 


meoe= transform of the driving function, f(t). 
Because the sound field is independent of the azimuthal angle, this 
transform can be formulated in the more manageable cylindrical co-ordi- 


nate system. In the new system the wave equation appears as: 


92P . 1 OP . 92P ae 6(Z-Zo) 
are a 5A a Bae + k*P = -F(w) Dar dr CaS 





where: 
u) 
kK = wave number, ; 
Cz) 
zZ = depth, and 
P = P(r,z,w), a function of horizontal range, depth and frequency. 


The Hankel transform allows one to compute the changes of pressure 
with a change of depth, independent of the range. It requires consider- 


ation of the vertical component of the wave number, kz: 
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ae kod = -s— 8(2-Z5) 2.4 


where: 


DB = B(E,rw) a function of horizontal wave number, horizontal range and 


frequency, 
ke = k?2 - &? 
—€ = horizontal wave number, and 
K = Wave number, AG 


This equation is then solved by variation of parameters. Let 


Bb = wi(zju@) * wae 2.5 
where: 
Zz 
= F Vil ZOU Zee2in) Lz oie 
Wi(z)=-— —— —— Sa 
, OT ee 
0 
h 
F Utz e(zazndz ; 
W2(Z) = 5— ; = iene 5b 
Way = U = a = (the Wronskian), Lane 


and both u and v are solutions of the homogeneous form of equation 2.4. 
Now if v(z) satisfies the upper boundary condition and u(z) satisfies 


the lower boundary condition, then: 
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Wy = OT Wuv 
O ZnS 2.6a 
( Zo 
and 
pF U(Z) OAT 
ieee TN 
° ow Nap 
0 ZL 


Hips means that the transformed pressure can be written in terms of depth 


dependent u(z) and v(z) as: 


Ge Zig) ca) 


Ee oe ie MO ueZe< ZG Za 
x Sey) (Ze) E 
Dp = oe aie am Ome tin e240 


A convenient way of representing this is 


~ . F Uulzy)v(zg) 
oe ae Zo 


mere: 
meececpresents the larger of z and z,, and 
mapeeprescnts the smaller of 2 and Z,. 


To obtain the actual pressure one must next take the inverse Hankel 


transform of zero order: 


ey 


BJA Er)edg 2 


Oo 

1) 
—_-_----— 
© 8 


2 ~ | meauaan dCem ede 2.10 
Use is made of the fact that [Ref. 1] 
Jee) = ee | H6 Cer) + ts? cer) 2 ell 
and 
Hoe Cerevin = -' (er) 2.12 


It is then possible to express the field integral (Equation 2.10) in 


terms of a Hankel funetion of the first Kind there 


ale \v(ze) 
a iF De! Hg | Cer) ede eae 
— Wuv 


The Hankel function can be approximated by 


Ws Cer) yf =2 ai era nee ar 5 
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Restricting this expansion to the first term makes the integral much 
more manageable but restricts the minimum horizontal range so that Er 


cE > |: 


Lye  ulzy)v(ze) j 
Pa y= e a a VE 2,15 


ae uv 


This field integral represents the pressure, within the limits 
imposed by truncating the expansion of the Hankel function. This integ- 
ral forms the basis for the derivation of formulae describing pressure 
meme CWO different numerical approaches. At this point, the approach to 
the normal mode solution, which will be used in the program EXACT, takes 
a different tack from the approach which will be used in the FFP solu- 


tion and program. 


B. NORMAL MODE METHOD FOR PRESSURE DETERMINATION 

In this paper the normal mode method refers to the process whereby 
approximations for the eigenvalue of a mode are successively refined by 
numerically integrating the square of a depth function across the depth. 
A correction term is developed from this integral and is used to refine 
the approximation until the correction is negligible. 

beuethe guess at tne mormal mode horizontal wave number, €,, be €. 


For §&, v(z) satisfies 


VT ks yo = 0 Ze oa 


i 


where: 
FOr Engen 2) ode lomlies 


Tae kee 2.16b 


where 
kgn = k* -En. 
Multiplying 2.16a by v, and 2.166 by v and subtracting yields 
Vav" > VV" + Wigs a Ke ee 2aly 


Then 


= (VaVi=Vv_') + Vig eat ee sO orn 


= (Va =a) =e eee vere 2.0 


Integrating from depth 0 to depth h: 


h nh 
(Vivi 7 Vvy") , = (&* -— €4) | VVy dz 2.19 
0 


or 
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h 


te Capea ev, (vi) = Wy00)v"(0) * v)(0)v(0) = (E?-E%) | Vvyndz. 2.20 
QO 


The following argument holds for any set of permissible boundary 
Pemearcions but, in order to clarify the discussion, the specific case of 
pressure release at z = 0 and rigid at z = h will be selected. The 
function v(z) is not a mode, but it can be forced to satisfy one of the 
poundary conditions without losing generality. if < Viz taseaer ined. vo 
Satisfy the boundary condition at z = 0, then, for any permissible boun- 
= evecondition, the last pair of terms on the left hand side of equation 


Bee0 1S identically zero. Therefore, 


h 


Ce emery (aati CE * 65) | VVp, dz. eae 
0 


Poreciic Simplest case considered in this paper, the bottom boundary con~ 


dition is rigid so that v,'(h) = 0. Therefore 


h 
Vyn(h)v'(h) = (&* - ES) | VVn dz eee 
0 
piideecince v= Vv, for € close to &,,, 
way Oh 
be? weheg = 


Zi 


Vic) 


fvdz vee 


i a 

Besides proving useful in determining SuccesSively better approximations 

of €), this same integral forms part of tne expression fom Uliemeon nome 

tion of a mode, since it 1S proportional to the derivative of thewiren- 
skian with respect to wave number. 

The Cauchy integral enables one to write equation 2.13 as Cheweaaiias 

the residues of all the poles of the integrand. Each pole corresponds 


£O a normal mode. 


(1) 
UnYyHo CéEnr)é 
pe a 2.21 
aT 
a; aw 
dé En 
where: 
WCEn) = 0 
To derive an expression for (dW/d&) multiply both sides by 
Un (hn) Ur! 
v,th) (8 also equals wnt 
h 
Un(n)v'(h) - ug(h)v(h) = =e 60)9 SZ Ques 
0 


The left-hand side of this equation gives the Wronskian for u and v, 


and, aS & > E,, WCE,) > ©, aS 10 should: 
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AW 


coli Sa) sl oe 
dé an At>0 Ag 2.26 
h 
WCE)-W( Ep) 
= = = B(E+En) wed, Py ik 
& Sn 
0 
Since —€ ~ En, and by properly constructing the program, the constant 8 
can be forced to be equal to 1, 
h 
dW 
a = 2 v*dz 2.26 
ae le, | 
0 
and equation 2.24 can now be written: 
ify ran (1) 
Past 2 SE ved ye Charen encod 


Using the exponential for the Hankel function, as in eqn. 2.13: 








_ iF 2 UnYn Gepr -i n/4 
Poa 2 Tenr Ivedz® © ae 
Se 
as Peek ees 74 ~ 
Yo 7 TE? [v*dz 


Zo 


In reality, this only represents the pressure in a system for which 
there are no boundary losses: i.e., for self-adjoint boundary conditions. 
If real losses, including scattering, were to be taken into account, 


nh 


the term v*7dz must be replaced by [Ref. 3]: 
Q 


h a 
di du 
v7dz + Hv2(0)($2) = ven) $2) eee 
| dé dé 
where: 
oe dv/ dz diay dz, 
Vv - u 


Ignoring boundary losses, equation 2.31 represents the total pressure 
at a point which is at depth z and horizontal distance r from a point cw 
source at depth z, This equation is the basis of the computer program 


EXACT, which is discussed in detail in Cnapter III. 


C. DIRECT EVALUATION OF THE FIELD INTEGRAL 

The second method of determining pressure at a point requires calcu- 
lating the field integral itself atte then numerically integrating. Equa- 
tion 2.15 represents the pressure field in terms of the Wronskian and 
the source and receiver depth functions u(z) and v(z). The integral is 
solved by replacing it with a summation across an interval of wave 
numbers. Pressure at a point is found to be the sum of the pressure 


contributions from each horizontal wave number: 
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én 
_ F 2 -in/t So UlZ>)v(zg) image 
Fe Te a e Pe ye e MAE AE Zs 


ae 


u(Z>), v(z¢), and Wy are calculated from the known boundary conditions 
at the surface and at the bottom. A&E must be small enough to ensure 
mecuracy. 

It is impractical to pursue this theme without introducing absorp- 
tion. Evaluation of equation 2.34 when mA&— = E, produces Wy, = 0 in the 
denominator, leading to an infinite solution. The magnitude of P, then 
depends on how near mA& comes to E,. Introduction of absorption as an 
imaginary component of horizontal wave number displaces the poles from 
the real axis and ensures that the Wronskian does not go to zero along 
the path of integration. This has the effect of dispersing the energy 
out of the theoretically limitless pressure function that results when 
the Wronskian goes to zero at the exact wave number corresponding to the 
normal modes. 

The terms within the summation sign of equation 2.33 have the same 
form as the discrete Fourier transform, normally considered for time, t, 


and frequency, w: 


N~-1 
G(nAw) senate = g(mAt). 2.34 


n=0 
By simply renaming the variables, the transform becomes: 


Z5 


N-1 


tes G(n AE) a oe = g(mAr). 2.50 
=(Q 


a 


This is the basis for the tecnnique known as the Fast Field Program 
[Ref .4]. By comparing equation 2.34 to equation 2.35, the following can 


be seen to be true: 
G(nAc) See 2.36 
W *) E=nAgE 


The pressure at range R can then be written as 


Ff 2. vin’ 
Pp = Paap = GeV oa -” Aeirenrent aNeT 


The capability of this Fast Field Program will be limited by the number 
of range (and wave number) increments. The wave number sampling incre- 
ment, A&, the range resolution, Ar, and the number of points in the FFT 


can be seen to be related by: 


-te Er = (nA&)(mAr) 2,38 


OF 


Teen Ne Ars ; 25a 
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morenm results from comparing equations 2.34 and 2.35. The number of 
Sampling points, N, places practical limits on either the maximum range 


or the range resolution. 


D. ABSORPTION AND BOUNDARY EFFECTS 

Sound energy is "lost" by absorption in the water and transmission 
Maceo the bottom. In both cases a primary variable is frequency but the 
angle of travel of the mode and the type of bottom must be taken into 
aceount. This paper limits discussion to a non-scattering surface, water 
volume and bottom. It also neglects energy that is transmitted into the 
boewom and then returned to aie water after refraction within the 
bottom. 

Sound absorption in the water at low frequencies (under 1 KHz) is 
due primarily to relaxation of the boric ion. Standard texts (Refs. 5,6] 


define absorption in this frequency range as: 


a = ae 2.40 
Where: 
a = absorptive loss in dB/m, and 
F = frequency in kilohertz. 


This loss ean be converted to the fractional loss, a, for a given dis- 


tance in nepers/meter: 


. 


Loss 


iT} 


20 logig ear dB 


ar = 20 logj9 e@ 


eOar logig e 


on 


Say are 


il 


f 
il 


8.7a dB/meter, or a = oa nepers/meter Cone 


This loss rate can be introduced into the equations derived earlier 


by adding absorption as an imaginary component of the wave number: 
Ka = K a la 2gie2 


The exponential form provides the absorptive loss factor: 


ai kok _ ai (k+ia)R 
TKR ek 
2 S ; 


2.43 


However, all of the expressions for pressure are given in terms of the 


horizontal range, r, not actual “path Miener a nee olnee 


xX hon 


a 
ae the loss 


factor can be written as: 


loss = e. ee 2.44 


Bottom interactions cause losses and phase changes which must be 
taken into account. In their discussion of propagation loss using normal 
modes ands an impedance boundary condition, Koch and Lindberg [Ref. 7] 
begin their discussion by defining the relationship between pressure and 


its derivative with respect to depth: 
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where 


f= tluid pressure in the water at the boundary, 


< = derivative of pressure with respect to depth, 


kK, = vertical wave number, and 
R = complex relfection (Rayleigh) coefficient. 
The Rayleigh reflection coefficient for the boundary between two 


homogeneous fluids is given [Ref. 5, Equation 6.30] by: 


PHoh o sin8p 











Pe ee one 
Pplh Sin@p, 
ee lh 
OP wly sine, 


where 
Cy - sound speed in the bottom at the interface, 
me eoound speed in the water at the interface, 
OP, = density of water, 
Ph = density of the botton, 
GO, = angle, measured from the horizontal at which the wave strikes 
une 


bottom, and 
6 = angle, measured from the horizontal, at which the wave is 
transmitted into the bottom. 


Substituting this value of R into Equation 2.45 and simplifying: 
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Oe sind 
- = j 7, eet s * D P, Ze Ane 
Z PDCph Sindy 


But 
Sin8@p = /1-C0S?6p 2.46 
Snell's Law relates the incident and transmitted angles by: 
oo 
cos 8h = = cosy ena 
Cw 
50, 


Ce 
SIN Opes ay |) ia - cos*6y . . 2S 
W 


The critical angle, 6,., is defined as that grazing angle for which 


the transmitted angle is zero. This means that: 





Sew, geos) and =). 2.51 


where 8, 1s the critical grazing angle, measured from the horizongeue 


Therefore, from Equation 2.50, 


cos*6,, 
sin6, = 1-7 2.52 
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by cef inition, 


cos 6, = 2 na 
and 


z Kg 
sin Ow = kK ° 2254 


Peieotloicing Equations 2.52 and 2.54 into Equation 2.47: 


dP ., Pwlw cos *6,, 
dz ~ \* PHCh y COS “6, 5 ey 


Im terms of wave numbers, Equation 2.53 can be used to give: 








dP ., Pwlw er 
dz ~ ** PHeD eee OS = 6. aoe 
or 
Pye 2 
ey 2.57 
dz PbCh COS “85 


Equation 2.56 shows that, for any grazing angle less than critical, 
the pressure derivative will be real whereas, for grazing angles greater 
than critical, the derivative will be purely imaginary. If the pressure 


were complex, the derivative would be complex and complicated, too 


syle 


complicated to solve in all cases with the basic program design. Therefore, 
only the real component of pressure was used in. the determination of the 
bomidary cCondlvions. 

A practical discussion of the physical Significance of the reflec 
coefficient and the relationship amongst it, the grazing angle and criti- 


cal angle is included in Chapter IV. 


a2 


Ti PROGRAMS 


Implementation of theory into the computer programs required two 
distinctly different methods. The program EXACT made accurate estimates 
of the normal mode horizontal wave numbers in order to compute the 
pressure as the sum of the pressure components from each of the modes. 
Absolute values of pressure are of no concern in either EXACT or FFP 
because the purpose is to compute propagation loss. It was assumed that 
forcing function, F, had a value of unity. The calculated pressure can 
then be compared to the pressure at one meter, and from that ratio, the 
propagation loss is calculable using decibel ratios. The program FFP did 
not compute horizontal wave numbers of modes; it assumed a large number 
of wave numbers equally spaced, calculated the pressure for each, and 
Summed the individual contributions to give a final pressure. 

Both programs were written using complex numbers and double preci- 
sion (64 bit words), except for subroutine ‘'eigens' which was restricted 
to single precision because it relied oem lMsoberouvine EOQRTIS which was 
available only in Single precision. 

Subroutines 'modes1' and 'modes2' provide many common functions for 
Beth programs EXACT and FFP. They calculate the transmitter and 
receiver modal pressure functions by stepping upward or downward from 
one of the boundaries, computing the pressure and its derivative with 
respect to depth by means of a fourth order Runge-Kutta technique 


Ger. 7]. 
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A. NORMAL MODE PROGRAM 

This program is attached as Annex B. EXACT computes the aporoxaae. 
modal horizontal wave numbers by means of the subroutine ‘'eigens'. tote 
then refines the horizontal wave number estimate to the desired accuracy 
using subroutine 'modes!i' or 'modes2' and computes the pressure faepors 
for each mode. Its aim is to compute the total pressure and, hence, the 
propagation loss at each range of interest: 

Subroutine 'eigens' uses a matrix method [Ref.8] to approximate 
ali the real modal horizontal wave numbers as well as the first evanes- 
cent mode. At the ranges of concern, it is unlikely that more than the 
first evanescent mode would be significant. The matrix of wave number 
estimates is found in a relatively short time and, although nev guar- 
anteed to be accurate, it is complete; no modes are skipped. 

The next step involves finding the exact values of the horizontal 
wave number for each mode. The procedure utilizes the Fact: thav iam 
each mode, both boundary conditions for a depth function (u or v from 
Chapter II) will be met. By setting the depth function or its derivative 
with respect to depth to a known boundary condition at the top or the 
bottom, it is possible to find a vertical wave number that will result 
in the boundary conditions being met at the the other boundary. The 
boundary with higher sound velocity will experience an exponential decay 
in the depth function. By starting at this = boundary, its a relative, 


simple matter to solve for vertical wave numbers that result in 
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Sinusoidal variations that lead to solvable boundary conditions as the 
other boundary is Eooneeeheds This avoids the possible problem that 
would exist if an attempt had been made to start at the boundary which 
mac the lower sound velocity. In that case it is possible that instead 
of solving for an exponentially decaying function, the function may be 
exponentially growing, and no solution is possible. Figure 3.1 illus- 
trates how a poor choice of starting boundary could lead to an unsolvable 
Situation. For each of the modes in the ‘main’ routine calls either 
'modesl' (when sound speed is greater at the top than at the bottom), or 
'modes2' otherwise. It calls with starting estimates of the modal 


horizontal wave number; the value from 'eigens' the first time and an 


improved value each subsequent time. 


'Modesi', using the best estimate for &,, computes the modal depth 
function at each increment of depth, starting at the surface where the 
depth function is zero and its derivative with respect to depth is set to 
1.0. The integral (Equation 2.31) is calculated at each interval and its 
mingle value at the bottom is used in the correction term. 

At the bottom, the depth function and derivative are compared to the 
boundary conditions expected for a mode. For a rigid bottom, a deriva- 
tive of zero is expected. For an impedance bottom, the expected deriva- 
tive will be defined by Equation 2.56. The error in the derivative, tog- 
ether with the depth function and the integral is used to calculate a 
correction for the modal horizontal wave Amise ee 

If the surface sound speed is greater at the bottom than at the top, 
Subroutine 'modes2' is called. It performs the same function as 'modes1' 


but starts at the bottom, where the depth function and derivative are 
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Set to meet the boundary conditions for an impedance bottom (Equation 
2.56). A correction is made after the subroutine has Breaced to the 
Surface and comparison has been made with modal boundary condition 
expectations; in this case depth pressure function equal to zero. 

"Modes1' or 'modes2' will be repeatedly called until £&, has been 
satisfactorily estimated. The ‘main' routine then computes the ‘incom- 
plete’ modal pressure from the modal depth functions, the integral, and 
bee horizontal wave number, —,, (Equation 2.31) for each mode, storing 
them for later use. It is termed ‘incomplete’ because there is no range 
dependency at this point. 

After the last mode has been processed by 'modest' or '‘'modes2', 
range is introduced and partial pressures for each mode are calculated, 
keeping an updated sum (complex) of all the modal pressure factors. 
These pressure factors are then converted to propagation loss. It is a 
Simple matter, once all the eigenvalues have been found, to compute the 
propagation loss for a series of ranges, aS can be seen in Appendix 


pee COgram EXACT. 


Be FAST FIELD PROGRAM 

This program was written with the intent of completing a Fast 
Fourier Program, FFP, which would utilize a Fast Fourier Transform, FFT, 
Subroutine. However, the program was only completed to the point where 
pressure factors were computed for individual ranges (as in EXACT), 
rather for a large number of ranges (as is the intent of an FFP). The 


computer program FFP, although incomplete, is included as Appendix C. 


on 


It was mentioned previously that EXACT and FFP are provided many 
common functions by 'modesi' and 'modes2'. There are also some basic 
differences in the way the subroutines are used by the two programs. In 
rFT, Since no exact wave numbers are to be calculated, there is no 
requirement for subroutine ‘eigens'. For the same reason, there is no 
requirement to find corrections to the wave numbers, so '‘modes1' and 
'modes2' are used differently. In fact, both Subroutines are required by 
cen | 

Using the incremental wave number, ndAé&, Subroutine 'modesi' starts 
at the surface with the same boundary conditions as in EXACT and steps 
down, computing the pressure function and the derivative, through the 
depth of the upper of the transmitter or the receiver, and stops at the 
lower of the two. Values are required at the lower level for later use 
with the values from 'modes2' to calculate the Wronskian. Then 'modes2' 
Starts at the bottom with the same boundary conditions as EXACT and 
Steps upward until the pressure function and derivative have been calcu- 
lated for the lower level. 

The ‘main’ Zacene then uses the pressure functions and derivatives 
to calculate a partial pressure (Equation 2.33) for that wave number in- 
crement. Each partial pressure is summed (integrated) so that, when 
the last wave number has been processed, the result is the total pres- 
Sure. 

As in EXACT, FFP computes the propagation loss for eacn range. 
Implementation of the FFT would make it unnecessary to step througn 
ranges aS must be done with EXACT; propagation loss would be calculated 


for as many ranges as there are FFT points. 
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IV. TESTS RUN ON THE PROGRAM 


Because it was difficult to gauge the correctness of the program or, 
more to the point, the veracity of this author's ideas, it was considered 
Beueent to check the results at each stage of development. Although the 
process was awkward and time-consuming it proved necessary and resulted 
in some unexpected benefits; the results of some of the tests provided 
graphical insight into the phenomenon of sound propagation. By starting 
With isospeed conditions, verifications were simplified. SOLULIOWS a or, 
horizontal wave numbers which resulted from the computer programs were 
enecked against ‘correct' solutions from simple formulae for the non- 
absorption, rigid bottom system, as well as for the system that included 
absorption in the water. Once the impedance bottom was introduced, an 
analytic approach was required to decide if corrections were applied 
properly. A simple expedient to cheek on the program at any point in- 
volved plotting propagation loss against horizontal range for pairs of 
modes and observing the interference distance between nulls. This dis- 
tance, when compared to the theoretical distance would highlight an 
error if there were an anomaly. Another check involved the observation 
that, for a shallow channel such as that used in the model, losses are 
generally spherical at short ranges and cylindrical for far ranges. The 
“smoothed propagation loss curves could be expected to lie roughly along 
curves predicted on this basis. 

The first tests involved studying the basic building blocks of the 


program. As the data were read from input files, they were printed into 


ag 


a new file to ensure that correct information was being utilized. 
Another simple check involved printing a matrix of the wave numbers as 
they were computed. To illustrate the importance of such a basic check, 
it is worth remarking that the half-increment method required that the 
wave number be calculated for each half-increment of depth, including a 
depth of zero. Since the variation of Fortran used does not supporta 
matrix with a zeroth element, a printout proved very useful in visualiz-— 
ing the situation and pin-pointing a subtle error. In many cases it was 
only by printing out all of the variables after each step that errors 


could be identified and remedied. 


A. RIGID BOTTOM WITH NO ABSORPTION 

Subroutine 'eigens' was first checked for correctness by comparing 
its constant speed solution, for each mode, to a solution which was 
known to be accurate. This accurate solution was based on tne concept 
that, for a constant sound speed, with rigid bottom and no absorption, 
the pressure function would be sinusoidal, satisfying both boundary con- 
ditions. This means that the horizontal wave number for each mode (n 
always odd) can be defined in terms of the wave number, k, and the water 


depth, H: 


pee 7 MU 
Calas (sq) *- Ho] 


Because the subroutine could only be run in Single precision, high reso- 
lution was not expected. Table 4.1 shows that resolution increased lin- 


early with the number of depth increments. For example, for mode 2, the 
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TABLE 4.1. SUBROUTINE 'EIGENS' - COMPUTED SQUARES OF WAVE NUMBERS AND 
ERRORS RESULTING FROM VARIATION OF NUMBER OF DEPTH INCREMEN hos 


MCDE NUMBER 
i 2 3 + 
NUM OF VALUE VALUE VALUE VALUE 
INCRMT ERROR ERROR ERROR ERROR 
OS Onemtee47o2 OF 1959ee25s2 0.1525131932 0.0237859429 
Meme COOlOG21s  -c 0009Ces2Gs  -0.0032877883 -0.06145329¢66 
50 Gee cee eeOmmoss@iseee  § ©.15087442289 =0.0209414381 
MOMeCGG4 56745 ON 504525355 20)6015490440 -0.0167259096 


100 Cec s5 47361 Orbos. Sey Om eoOatOee7 §=0.0305268121 
memeOOG7Z 2775-0 0002279493 -0.0008256757 -0.0071405356 


200 O22143832949 Og 52023 OFte3 6225239 - -0.0342555159 
mero OOOtlee 7 oO =e, COClls6725 =-0,0004131089 -0.0033718318 


400 CE Leow 7 2 0.1951455441 CP lel 45202495 )-C . 0560238122 
mere 000057292 50-0 09000552758. -0.0002056200 -0.0016435355 


800 0.2143746840 Opie foi lidh i hess} Ore z2 273107 =0. 0368554579 
pOmovO@O22 7/60 s—-0,, 0000284275 -0.,0001033zZol1 -0.0008118898 


mows O.2143732463 OmlIS1O2Z59¢s GeeI2770720 =0.0372637914 
-0.0000014384 -0.0000142262 -0.0000516671 -0.0004035562 


moO 860. 2143725272 Or 509 Saez 1 Cees 22 oe -0'..03 74661565 
Orme OODOe7 193. C7 OC OCG Ia 7)6 = 0 200002583456 -0.0002011912 


mescT ©.2143718079 O-1SS0SS055 5 Orles2 25-049 -O.03766/3477 
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error using 25 increments was .000900, and the error was halved by in- 
creasing the number of increments to 50. This trend continued through 
the run that used 3200 increments and resulted in an error of 0.000007. 
it was decided that an error of 0.00002, whicn resulted from the use ver 
100 inerements, was acceptable. The value of each norizontal wave 
number determined in subroutine 'eigens' was intended only as a starting 
point for the more accurate subroutines 'modes1i' and 'modes2'. The 
accuracy, therefore, has significance only in that, if two modesyyag- 
closer together than single precision accuracy limitations, there would 
be a chance that a mode will be missed completely. 

Subroutine 'eigens' was also checked using sound speed oror ites tnae 
varied linearly, botn positively and negatively, with depth. Since no 
clear-cut means of assessing the error was available, no conclusions 
could be drawn, and results are published for only the positive gradient 
(Table 4.2a). In addition, by using a sound speed profile that produced 
two strong ducts, it was possible to show that 'eigens' was capable of 
distinguishing between two horizontal wave numbers which are very close 
together (Table 4.2b). The solutions are not precise but no modes are 
missed and the necessary information is provided to the subsequent sub- 
routines where the values are calculated to the required precision: 

Next, normalized modal pressures at a large number of depth incre- 
ments were calculated using 'modes2'. These data were combined in 
Figure 4.1 and show, for each mode, the relative amplitude of the modal 
pressure function. The profiles for all three real modes as well as the 
first evanescent mode are plotted together and care must be taken to 


note that they represent the modal Characlerisvics enev vie wacrual 
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pressure contributions. The figure shows that, for rieid bottom, no abaern— 
tion, iso-speed conditions, the profiles meet the boundary conditions at 

Che top, where pressure is zero, and at the bottom, where the derivative 

is zero. The same test procedure was repeated using subroutine 'modes1' 

and identical results were obtained. 

To illustrate the relative degree to which each mode is stimulated, 
depending on the depth of receiver and transmitter, a trial run was made 
With a transmitter at depth 13 meters and receiver at various depths. 
Figure 4.2 depicts the relative strength of each mode at each depth. 
The second mode is most strongly stimulated and mode 1 is least 
strongly stimulated, as could be anticipated from Figure 4.1. 

Routine FFP does not use exact normal mode horizontal wave numbers, 
but boundary conditions must remain consistent and sinusoidal variations 
still occur at a rate defined by the vertical wave number, kz, whieh, in 
turn, is a function of the horizontal wave number, &. The results of 
subroutine 'modes!' were compared to the values of sin(kzzy) and the 
results of 'modes2' were compared to the values of cos(k,(H-z<), whercu 
is the water depth and zs and 22 \anemei- mcereicmme 6 ec receiver and 
Cransmitter. The comparisons showed that the subroutines worked 
properly. 

At this point, the solution using calculated horizontal wave numbers 
in EXACT was examined for accuracy. It was decided that the accuracy 
criterion would be based upon the worst case which allowed for all 
errors to be cumulative. It can be seen in Equation 2.31 that an error 
in E, will exhibit itself directly in the Square reousand in tne exeon-ms 


tial. A small error causes only a small error in the square root but it 
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is critical in the pnase component. It was categorically decided that 
the maximum allowable phase error for each mode would be 0.1 radians, 
and that the maximum range would be 100,000 meters. This limits tne 
maximum error in &), to @ 2torvall wave: numbers... 

In assessing possible real errors, a vector plot was made for each 
of several solutions. A constant adjustment was then introduced to eacn 
@r the horizontal wave numbers, giving them an artificial error. It can 
be seen in Figure 4.3 that for short ranges, a large error is admissible 
and, even at 50,000 meters, the limit imposed by the accuracy criterion 
is well within the allowable error. 

Vector plots were also used to illustrate how a small cnange of 
range can cause a dramatic change in pressure. Figure 4.4a shows the 
change of propagation loss from 830 to 840 meters. Most of the change 
is involved with phase and the amplitude change is small. This can be 
seen best in Figure 4.4b where the range change i8 only one meter. It 
should be noted too, that these errors have been made cumulative 
whereas it is unlikely that the errors would ever be that way. One 
would therefore expect the total error to be less thane enar 1llus- 
erated. 

A plot of loss with range for a rigid bottom, Figure 4.5, showed a 
general 20 log R increase in propagation loss close to the source and a 

. e 
change of 10 log R as the distance increased. Superimposed on this is a 
strong interference pattern from the three real modes. Sudden changes 
of pressure with changing range, as were noted in Figure 4.4, are seen as 
Piolo lLuGhualions Or propagation loss in this plot. A decisive check on 


Che general reliablity of the program to this point involved excluding 
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all but two real modes and plotting the resulting propagation loss 
curves. For each pair of modes, the interference distance can be shown 


bo be related to the horizontal wave numbers by: 


eee al Wee 
om Sn 





alee 
R is the horizontal distance between nulls, 
and €, and _ agesehe rédl horizontal wave numbers. 
Comparison was made between the calculated interference distance and 
the observed distance for all three interference combinations and proved 
to be correct in all cases. Figure 4.6 shows the result when the first 
and third modes were used. The cycle distance was calculated to be 96.5 


meters, which agrees with the distance from the plot. 


B. RIGID BOTTOM WITH ABSORPTION 

Absorption due to water (low frequency), although normally very 
small, was made artificially high to cheek on the mechanics of the 
program. With absorption set to an arbditrary value of 0.0005 
nepers/meter, it was noted that the imaginary component of the wave 
number increased with the increasing grazing angle associated with higher 
fmes: 0.000505 for mode 1, 0.000555 for mode 2, and 0.000727 for mode 

To ensure that output losses were consistent with those input, other 
Beeparacion loss curves were drawn uSing a variety of rates of absorp- 
tion. Figure 4.7 shows that the difference in loss at any given range is 
8.7 times the difference of absorption rates (as expected from Equation 


Bee) For example, a change from 0.0000001 nepers/meter to 0.0001 
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nepers/meter causes a change of propagation loss, at 10 kilometers, of 
paoyde. One would have expected a change of 1.0 nepers or 8.7 dB. 

A series of runs was made to illustrate how pressure varies with 
receiver depth at a set range for a variety of source depths. From tne 
grapns (Figure 4.8) it is important to note that for a source or a 
receiver at the surface, the pressure will be zero and the propagation 
loss infinite because none of the modes is stimulated. This does not 
change when an impedance bottom is introduced, but it would change if a 
real Surface were considered. 

Using these same depth profiles, it was Merete 65 further verify 
the program by checking for reciprocity. The propagation loss is, for 
example, the same for a combination of a source at 5 meters and receiver 
at 25 meters as for a source at 25 meters and receiver at 5 meters. 

The final set of tests for a rigid bottom involved the use of the 
FFP. At the point of testing, the Fourier transform (FFT) had not yet 
been introduced so pressure factors were being calculated for individual 
Prices. Figure 4.9 shows that, for a zero absorption loss, the only hor- 
izontal wave numbers that contribute to the total are at the exact mode 
values. Sieameemchneaniges Of the “forizontal) wave number cause Serious 
changes in the value close to the exact mode value. 

Introduction of absorption causes the energy to be dispersed across 
the spectrum. There are small contributions from across the wave 
number spectrum and flattening at the modes. Higher absorption causes 
Meeaver dispersion of the ‘pressure’. The integral, with absorption, 
can be seen to fluctuate in Figure 4.10, staying close to zero until the 


first mode is approached. (jemi e te emimeteases Gulickly, but not 
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Figure 4.9. Amplitude of Sound Presstire Factor for [n@ivagiuall 
Wave Numbers with No Absorption. 
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Figure 4.10. Integrated Pressure with Absorption Rate 
0.00005 Nepers/Meter at Range 10000 Meters, 
Tx/Rx at 10/37 Meters, 50 Meter Ocean, 50 Hz. 
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abruptly when a mode (5 eee A sharp decline is notable when the 
180 degree phase change occurs at the mode. The final value represents 
the final pressure at the receiver. 

Various checks were made to ensure proper functioning of the FFP 
program for a given range. By changing the number of wave number incre- 
ments, it was found that the maximum amplitude for each mode changed 
(Figure 4.11) but that the integrated pressure was identical when the 
number of increments was changed from 1024 to 2048. Pressures agreed 
With one another to the seventh decimal and fifth significant digit 
Given more time, it would have been interesting to see how few samples 
could be taken, for a single range, before significant errors were intro- 
duced. 

Figure 4.12 shows the effect of changing absorption rates. However, 
the change of pressure at a given range does not correlate with the 
Change of absorption rate. This indicated that the FFP routine was in- 
correct in some way and further checks were required. 

The integrated pressures from this program were then compared to 
the values obtained from the program EXACT for different absorption 
rates and ranges. As can be seen in Table 4.3, agreement was poor. 
Because of time constraints, it was decided to abandon the FFP and con- 
centrate on solutions using the EXACT method: 

NOTE: A great deal of time had been spent on FFP trying to find a solu- 
tion: £Cr Vamaxead ‘range, a prerequisite to implementation of the ‘ERP. 
subroutine. At the time that work with the FFT was first suspended, it 
gave inconsistent results that did not meet any of the testing criteria. 


A small change of increment size, from 2047 to 2048 increments, caused 
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Meee 4.3, CALCULATED PRESSURE FOR PROGRAM FFP AND EXACT. 
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Figure 4.12. Integrated Pressure for Absorption Rate 
a. 0.00001 Nepers/Meter (Lower) and 
b. 0.0000001 Nepers/Meter (Upper) Range 1000 Meters. 
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outlandish changes in the integrated pressure. The changes of pressure 
resulting from a changed absorption rate did not correspond, even rem- 
otely, with the cnanges expected. At that point, op ions were shifted 
and it was decided to concentrate on solving the ‘EXACT' program prob- 
lems at the expense of the 'rTFP'. At a later date, too late to be prac 
tically helpful, some of the FFP problems were rectified and the above 
mentioned inconsistencies were remedied. It’ was possible Ut0) geiiae 
propagation loss to absorption rate and it was shown that the number of 
Samples has a negligible effect on the final pressure. However, the 
discrepancy between the EXACT solution and the FFP solutions still exist 


and the FFP is incomplete. 


Oi 


eee PEDANCE BOTTOM 

Tae introduction of an impedance bottom to replace the perfectly 
reflecting rigid bottom made it possible to better predict the propaga- 
tion loss experienced in reality. This was done by defining tne bottom 
in terms of its density and the sound speed in the bottom at the inter- 
face. After continued failure of the program to solve for modes whose 
grazing angles were far below the critical angle and the eventual reso- 
lution of that problem, a series of tests was devised to try to optimize 
the program. These tests also provided a convenient means of illustrat- 
ing how various angles of incidence associated with the different modes 
resulted in the varying phase changes and amplitudes of reflected 
energy. 

As with the rigid bottom case, it was possible to plot the modal 
pressure as a function of depth, from the reflecting surface to the imp- 
eeance bottom. Figure 4.13 shows how the pressure reaches a maximum 
well above the bottom for mode 1. The grazing angle, in this case, Is 
less than the critical angle and illustrates the pnase advance at the 
bottom. The test was repeated for strong positive and negative gradi- 
ents bUcC Showed little change. 

The first set of tests to check on final accuracy involved using a 
critical angle of 0.3 radians and was intended to find the optimal 
number of increments. It had been determined earlier that accuracy to 
the ninth decimal, or better, was obtainable (iso-speed) for all modes 
using 200 depth incements/mode. At 100, the daecuracy dropped to the 


seventh decimal and at 50, answers were accurate to only the fourth 
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decimal. Table 4.4 shows horizontal wave numbers that were calculated 
Wegme 100, 200, 400 and 600 depth incr ements/mode Mie COmbpination Wilh 
@erereccion factors of 1, 2, 5 and 10. The correction factor was divided 
into tne mode wave number correction term of subroutines 'modes1' and 
miedeac'to limit the size correction and thus prevent instability 
These correction factors were found necessary for modes whose grazing 
angles were less than the critical angle (in this case only mode 1). 
Without perc OLGCe slemeracvor, one program would fail to converge to an 
a@eeceprable solution but would become divergent. 

Peeiajor problem in deciding the best choice of variables came about 
because there was no correct answer against which to appraise the 
values. It was reasoned that 800 increments per mode and a correction 
factor of 10 would give the most accurate answer and the problem nar- 
rowed down to finding a combination that gave an acceptable accuracy 
with a reasonable number of iterations. Having decided that 200 incre- 
ments/mode was sufficient and necessary for the rigid bottom case, it 
became a question of how large the correction factor could be. ihe 
Smeree Was further narrowed to a correction maeceer Of 2 with an error of 
9.8 x 10°® and 40 iterations required and a factor of 5 with an error of 
3 x 10°° and 89 iterations. Because the error was borderline for the 
@emreection factor 2, the final choice was 200 intervals/mode with a cor- 
meeeron factor of 5. 

Using the inputs as decided above, an effort was made to correlate 
actual accuracy with that specified as a requirement. Results for 


trials on all modes using accuracy criteria from 1075 to 107*? follow in 
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ee 


TABLE 4.4. PROGRAM 'EXACT' - REAL AND IMAGINARY COMPONENTS OF € 
FOR MODES 1 TO 4 RESULTING FROM VARIATION OF NUMBER 
OF DEPTH INCREMENTS AND “CORRECTION” —EACIO@rs. 

COR INC ITER 2 2 3 “ 
1 100 $0 0.212512834 0.195507949 ©. 150035436 C.002867414 
0.000450164  -0.C0333983i  -0.910172100 0.0377799m7 
2 100 38 0.212512965 0.19550794¢ 0. 150035236 9.c02sé7424 
0.000450172 -0.003339831 -0.010172100 0.037779947 
5 100 289 O221251905.- 0.29550794¢ ©. 150035436 0.0025574i4 
0.000¢50172 -0.003339831  -0:910172100 C.037779947 
10 100 285 O.212513060 0.195507949 0.150035436 0.002867424 
0.000450170 -0.003339831  -0.¢10172100 C.037779947 
1 200 74 0.2125128s5 0.195507947 0. 150035443 0.002967455 
0.000450155 -0.003339830 -0.c10172094 C.027779560 
2 20Guures 0.21251296a 0. 195507887 0.250035443 c.co2semess 
0.0C0450171 -0.003338830 -0.010172094 ¢.037779860 
5 200 §3 O.212513035 ©.195507947 Q. 150035443 0 .COmeemee 
0.000C+50172 -0.00322983G -9.01017209+ 0.723777 2360 
10 200 195 0.212513063 0.19£507947 0.150035443 0.COS6gaes 
0.0C0450170 -0.06333¢830 -0.010172094 C.037779680 
1 400 67 0.2125128s 0.195507946 0.150035443 0.002867+37 
0.0C0450165 -0.0033393830 +-0.016172093 ©.037779642 
2 400 40 O.212512953 0.195507946 0.156035433 C.CO2857 2am 
0.0C045017. -0.003339830 -9.610172¢93 C.C3777se= 
5 400 89 0.21251303 0.195597946 @2150055443 C.CO2S6 eam 
0.600650172  -C.003339830 -0.010172093 OC. C2777 eee 
10 400 136 O12125120s 0.195507946 0. 150035443 C.002557aam 
0.0C0450170 -0.003339330  -0.010172¢093 C.C3777Ga=. 
1 800 72 O.212512685 0.195507946 G2156035 283 0.00286 725m 
0.000450165  -0.0063339830  -0.020172692 0.0377 7a 
2 900 39 O.21251296¢ 0.195507946 0.150035+43 C,cO28a7eem 
0.00045017: -0.0033329839  -0.010172¢92 0.037779442 
5 800 30 0.212513039 0.195507946 0.150035443 0.C02867237 
0.000450171 ~-0.003339830 -0.010172092 0.0377 74ae 
10 800 187 0.212513066 0.195507946 0.150035443 0.002857437 
0.000450170 -0.003339830 -0.010172092 0.037779441 
EXACT 0.214371820 0.195088818 0.149226333 0.002867437 
(RIGID BOTT) 0.000505340 0.000555289 0.000725949 0.0377794=0 


66 


Table 4.5. Assuming that the answers obtained using 10 *? had the least 
error and could then serve as a standard, errors were calculated as the 
difference between the standard and each answer. It was found, in all 
Cases, that the actual errors were slightly less than those specified by 
the accuracy criteria. This was to be expected since corrections were 
made to make the error less than the designated amount. 

The final check involved the problems associated with low modes 
whose grazing Beis were less than the critical angle. A series of 
calculations was made to find the number of iterations required for 
different critical angles (different sound speed ratios) using a 100 
meter deep, isovelocity ocean that supported seven real modes. It was 
Seoeericended to uncover further problems associated with large critical 
angles. At the same time, the squared values of horizontal wave numbers 
imermeune seven modes were calculated for each critical angle and are 
Seew in Table 4.6. For angles close to the grazing angle for the first 
mode, 0.1 radians, the program failed using both subroutines 'modesi' and 
'modes2' during these runs. This was unusual because previous data had 
peeomeolrcered for almost identical conditions. | 

The fraction of energy reflected and the phase change represented by 
the Rayleigh coefficient, were calculated for each mode for each criti- 
cal angle and are shown in Figure 4.14 and 4.15 (upper). The graphs show 
that the loss per bounce is zero when the grazing angle is less than the 
Critical angle. Mode 1 has tne lowest grazing angle and it can be seen 
that the critical angle must be very low if the low modes are to suffer 


any attenuation at all. Higher modes suffer higher losses. To show the 
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Figure 4.14. Phase Change/Bounce for 3 Real Modes as a Function 
Of Critical Angle. 
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extreme modal dependence of attenuation, Figure 4.15 (lower) displays 
the loss per meter for each mode for a variety of cirtical angles. 
These curves illustrate that the increased grazing angles associated 
Paeeoeeceniec higher modes result, not only in higher loss per bounce, but 
also in an increased number of bounces. For a small critical angle, say 
0.1 radians, the loss per bounce is approximately three times as great 
for mode 2 as for mode 1, but the loss per meter is about ten times as 
great (Figure 4.15). 

It was found that a change of critical angle (change of the sound 
speed ratio) changes the loss greatly. Figure 4.16 demonstrates the 
change of pressure due to mode 1 when the critical eee is changed from 
above to below the grazing angle. Changing the impedance ratio by ten 
per cent caused a corresponding increase in the loss but did not change 
Maemeoinit at which the critical angle became effective. 

Mmemoulldy the bottom loss effects for grazing angles greater Chan 
Critical, different runs were made for critical angles that allowed a 
different number of modes to propagate in the low loss manner (grazing 
angle less than critical angle). Figures 4.17 and 4.18 show how much 
each mode contributes toward the total pressure, depending primarily 
Mpen whether its grazing angle is above or below the critical angle. 
They show that for each mode, when the grazing angle is less than the 
critical angle, the propagation loss for that mode is identical to that 
for the rigid bottom case. 

As a last demonstration of the effects of the impedance bottom on 


propagation, data were compared for the isovelocity case and the strong 
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Figure 4.15. Loss Bounce (Upper) and Loss/Meter (Lower) 
Due to Bottom Interaction. 
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positive gradient case when the critical angle was 0.1 radians. It can 
be seen in Figure 4.19 that the upward refraction casued by the gradient 


lessens loss due to bottom interactions. 
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V. CONCLUSION 


The program EXACT, which computes propaation loss by first solving 
for individual mode wave numbers, was found to give realistic values. 
The Fast Field Program, which was meant to solve for loss using an FFT 
was not completed because of complications that proved unsolvable in the 
limited time available. 

One of the problems with the solution of the program EXACT, which 
would also be a problem with FFP, is the restriction placed on thea. 
dance bottom. All of the sound transmitted into tne bottom is eeeneae 
dered lost to the bottom. No allowance is made for this energy to be 
refracted upward and re-transmitted into the water. This might not be 
entirely unrealistic in some cases, because bottom absorption is often 
very high. However, this simplification limits tne practicality of tnis 
program and is considered by the author to be its main weakness. 

Another point of consideration is tne number of modes involved. The 
tests run usually involved only three real modes, two of which are lost 
to the bottom after a short range when the critical angle is 0.1 radians 
(sound speed ratio is 1.005:1). An increase of depth or frequency will 
mean an increased number of modes. There might be 4, 400, or 40,000 
modes present but a large portion of them may suffer a large bottom 
loss. This also explains why trapped modes are so significant and why 


it is often sufficient to consider only a few modes in most analyses. 
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One of the intended purposes of the thesis was to compare the com- 
puter processing time for each program to solve for conditions where a 
large number of modes were present. It was reasoned that an increase 
of frequency and/or depth would require a proportionate increase in tne 
mumoer of depth samples used by tne subroutines that calculated pres- 
sure as a function of depth for each horizontal wave number. This would 
mean that if the frequency were increased from 50 to 100 Hertz and the 
deoth from 50 to 500 meters, there ce be a 20 fold inerease in the 
processor time used to do these calculations in each of EXACT and FFP. 
fewever, there would be a further 20 fold increase in time for EXACT 
because it would be required to calculate for 20 times as many modes. 
The FFP would not require any similar increase and would therefore nave 
a decided advantage when a large number of modes were involved. [It is 
possible that for cases with a small number of modes, PRACT would be 
faster, but the information was not available for comparison. 

& Simplification in EXACT involves the accuracy criterion of the hor- 
Peelal wave Number. It was decided that errors due to accuracy limita- 
tions would be acceptable if the wave number were resolved to less than 
nor’, However, the subroutines ‘eigens', 'modesi', and 'modes2' return 
Mieesqliare Of the wave number. For the calculated & to remain within 


the limit, 
GA upemnor | < (Ey 107°) * 


SO 


tec | sGree aims En. 
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For the conditions considered, the correct accuracy criteria should have 
been marginally more restrictive for all three modes considered. Had 
any very small wave numbers resulted, the allowable accuracy criteria 
should have been much more restrictive. A more complete program would 
calculate eacn wave number to its own allowable accuracy limit using 
formula 5.1 

While it is recognized that program results would in no way alter 
the facts of nature, the sorreaseneeae between observed phenomena and 
predicted results did add to the credibility of the theory and program 
outlined in thls thesis Many of the results, as illustrated 1ng7— 
graphs, would prove useful as teaching aids to elucidate the principles 
of sound propagation in the ocean. Without reference to tne development 
of normal mode theory, it is possible to illustrate the consequences of 
Changing depth, range, frequency, water depth, and aeune speed profiles. 
For example, the program could be used to show-thne effects of surface 
duct on long range propagation. 

Finally, it is realized that the programming methods used are far 
less efficient than many readily available programs. MThe intention of 
this thesis was to develop an individual Normal Mode program from a the- 
oretical basis without direct reference to other, more sophisticated pro- 
grams. A more complete study would compare program results to actual 
observed data as well as to results of other, more authoritative 


predictions. 


80 


APPENDIX A 


MATRIX METHOD FOR FINDING MODAL HORIZONTAL WAVE NUMBERS 


The matrix method (Ref.8) employed in the program EXACT yields a 
Semplete set of eigenvalues in a Single pass utilizing IMSL routine 
EQRTIS. It is convenient, fast and has the advantage that all modes are 
measur edliy determined, no matter how close together they are. Accuracy 


is not attempted; that is assured with the Runge-Kutta routines that are 


used subsequently. 


The Matrix solution is a means of solving for the eigenvalues of the 


one dimensional differential equation: 


where: 


ols the eigenfunction, 
k, is the vertical wave number V(k*&°*), 


— is the horizontal wave number, and 
Z is the depth. 
It solves the equation by finding the eigenvalues of the matrix M in 


mae following Matrix equation: 


my = -(Az)* Ey ne 
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This equation results from writing the differential equation as a finite 


difference equation at each of N deptn Az apart. a column matrix with 


N elements,.eachn element definec by: 


Wn = W(ndz) Gee 1 fee A.3 


and M is a Square symmetric tridiagonal matrix which allows a computetion- 


ally efficient means of finding the eigenvalues. 


? 


Ds = 0 = a = 0 
=| OF = = = = = 0 
0 7 oe = : : 
Mos - = : = = = = ALY 
0 é = : = = Diy 1 = 
0 = a 7 = 0 = Dy 
J 
where: 
Da = 2 > | (hzjek Giz) Nee iy: 2 ee Ae5 
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APPENDIX B 


COMPUTER Sr ROGRAM EXACT 


See oe See : Deen il 

Ce ee RCGre eee abo UF TO 30 PAIRS OF SOUND VELOCITY INFORMATION ~°* 
Ce Piso ten menos tED VALUES OF WAVENUMBER FOR BEACH INCREMENTAL ° 
Ge Dette Peenuolun hs Lik NUSBER OF REAL ZIGEMVALUES AS WELL AS i 
(Ce MAX WAVENUMBER AND THE DEPTH INCREMENT Fs 
C= PESCR Ete, [5 nASoUMED-ICO Es CONSTANT AT ALL DEPTHS : 
Sai MON-COMPLEX EIGENVALUES ARE CALCULATED BY A MATRIX METHOD iN : 
ot Seee me eens AbD TAESE ARE USED AS A FIRST ESTiflaTE I¢ CLE » 
ee OF MODES SULROUTINES e 
oA Pro oURn= Ss cr > TRE WARONSKIAN ERC MODES ARE USED TO CALCULATE = 
Ce COMPLEX INDIVIDUAL MODAL PRESSURES WHICH ARE SUMMED FOR ‘ 
ee PARTICULAR RANGES a 
Ce BOLTOM ABSORPTION NOT ACCOUNTED FOR 4 
Ce i 


Sy a ne eae 


Se ooo oo ( LOOOl) CALe (a0) ARGUE, VAL, CORR, (NT,2DIES,FTS, 
ee Oe er Pe aoe oO yh lG (20) PRess( 50), PRESSR, IMAG, ALFA, AIG(50), 
Peat Ree, COMP 

Pe el oO) etme) pn ooo C130), 81G2(50),hMAG(50),NNAM, Ki, 

ieee WelOOOl Ra 10001), PRC SC), PiIlM(SO) 

ECGICAL. TCL 

TOP = .FXLSE. 
NUM=1 
Pig=j=03. 121592654 
CMAZ=O 

ChiITIN=200 


C2 


SUMP = 
ERRMAX 
CORN = 2 
Carrs 
RATIO 
ie =. Dit 
e 


G tl Il 
banal 
(a) 
— 
ead 
i, 
me 
~—— 
Q 
_— 
ss 
eG 
Ge 
~— 


Seren) er cee ee: 


QR RRR RR RRR EKER RRR ER RRR RRR RRR ER RRR RE RR ES Ree 
S 
IF (D(NUM) .LT. 0) GOTO 14 : 
KAY(NUM) = 2 * PI * FREQ / C(NUM) 
IF (C(NUM).GT.CMAX) THEN 
CMAX=C (NUM) 
END IF 
If (C(NUM).LT.CMIN) THEN 
MMAX=KAY (NUM) 
CIN = "ern 
END IF 
IF (NUM .EQ. 1) GOTO 13 


RATE (NUM-1) = (KAY(NUM) - KAY(NUM-1)) / (D(NUM)-D(NUM-1)) 
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1¢ 


es eee), coe 
PRINT 1039 D(NUi=1), C(NUH=1) RAY (CN) RATE ero 
NUM = NUM + 1 
COFG 12 
BOTTON=D(NUM-1) 
N = TE IXUCESBOTIO NR ERE tee = A220) 
WRITE(+4,6389) N 

6389 FORMAT (//' THERE ARE ',16,' REAL MCDES'/) 
M=1+N 
DELZ = BOTTOM/(INC*N) 


ee 


C(NUM-1) 
/CCNUM=1), KAY (NUMH-1) 
TOP = . TRUE. 


PRINT 103,D(NUM-1), 
WRITE (4,103)D(NUM-1) 
LE (C(NUM=1) LT ean 


IF (N.GT.S0O) THEM 
WRITE(4,100) N 
i100 FORMAT('TOO MANY MCDES FOR THIS FPROGRAM!! (REAL) = '.i4) 
COTeeadL 
ENDL 
c 
CTR RR EKKR RRR RTARTA TRE ER TR KK TR KE RR NRK eee ee YR eK ee ee 
CcC* 7 
Cw TRE NEXT SECTION CALCULATES THE SAnC i ie Sige en so 2On fee im 
oe ROUTINES EICENS, HODES] AND HODESZ 
oe THE REAL COMPONENTS FOR MODES! AND MOCDESZ ARE CALCULA 730 sae. > 
Ce INC*N + ] DEPTHS STARTING WiTH DEPTH 1 Bl [HE SURESCE Glo eee 
[a AL THE Bert : 
AT REGULAR DEPTH INCREMENTS A REAL VALUE CF K IS TAKEN FOR * 
ct USE IN SUBRCUTINE EIGENS THE NUMSER OF SAMPLES IS iCO : 
Cr 
CREVN TEN KKKKRK TORR TAR CRN RRR RR OE Se OR Re ne oe ae me Fe ae re 
c 
NOU = — 
NUM = 1 
Doers 1G = 1, 1Ne-tes 
De = (1G=1) 9 bekz 
KSsA(IG) = KAYQCNUM) + RATE (iG) (er =) 
CHFF RARER ARES RAAT AREER ESA RAR < 
C IT 1S POSSISLE TO iHSeRyY A LoOr QC reer ee 28 ee ee 
Gs THE DESTH Sar:PLeS OF K( ) A Test KRUPA = wee CC eee 
Cc VALUES OS PALNED FOR ican ae 
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18 


Se pL atte BD be tie ee a or 


le (INC. LE.200) JUMP = 2 
IF (MOD (IG-1,JUMP) .EQ. 0) THEM 
KOUNT = KOUNT + 1 
IE(IG.EQ.1) THEN 
KOUNT = 0 : 
GOTO 886 
ENDIF 
KAE(CKOUNT > = Kanga) 
END IF 
IF (DP.GE.D(NUM + 1)) NUM = NUM + 1 
CONTINUE 


CALL. EIGENS (BOTTOM, KCUNT MU PAE, ElG2y ality) 
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( 

a] 
“| 
. 


eh te) RON ce 
ee ee 


Ce Re ee ee ee Ree eee EK RAT KEK RAR NS RK © 


Gr ee ee ee ee Roy We Ea C..LCULATSD AS FUNCTION OF FREC i 


FREK = FREQ/1000 - 
FREEK = FREK**2 
Popeyes wOGCOL*EFREEK/( 1+FREER ) 
& 
Crewe eke kaw ee kb nn wk SRK KK ARERR ERT KN hake ee ke RE ET 
Gr i enero OCr NAIC Gail Eee USSD TO DETERMING THE PRESSURE 
ee ORNs TenGetICN AL ALL DEPTHS FROM TOP TO BOTTOM 


CER KRERRK RK KRAREKR RHR RRR ENTRAR RRR TRA RRR R RK RRR RK RR TA KK KT 


ALEAB = On 
ALEAL i W + ALPAB 
ALFA = CieE(0-¢, ALEAT) 

5500 DO-S9Se1k =41, Ilene 

RUS) a OACih) = -AEEn 

598 CONTINUE 
ARGUE=K(20)*72-(PI*=(I-.5)/BOTTOM)**2 
CALG( 1) = "@c0RT(ARGCUE) 


0 
Ee 
LX 


ec CALC IS USED TO EVALUATE THE CORRECINESS OF THE PROGRA! = 
Ce BY INPUTTING AM ISOVELOCITY PROFILE IT IS POSSIBLE TO CALCULATE + 
e- THE CORRECT HORIZONTAL WAVENUMBER FOR EACH MODE AND TEEN COMPARE > 
Cx THE RESULTS OF THE PROGRAM TO THIS 'CORRECT' ANSWER : 
Cee Re AR RR NO Ba i RRR eR He NR ee ee eo Ce a So oa 

VAL = EIG2(1) 

CORRN = 1: 
‘e 

Ge(eCPy THEN 


ee a ee Ke he eK eA KEARSE A TRE SS RT 


Cz Poe ool CeO weet EGP THAN AD Tae BOTT: 
or ee sPenOUt Meee se oe sb ou reenwisr MCDES2Z 1S USED 

Ce ene eee en ow Ceo onRO We tiPoe te CORRECTION MADE TO THE 
Cs SOUAR® OF THES HOs I ZO teyeavVeursGl 1S COMPARED TO A ERROR 

Ce rete Olle oe ee eee eee eee tilt tae SUBROUTILE IS CALLED 
CK AGAIN. (Rete eC ce Ree Pte lowteee Tell IRE LIMIT, THE PROGRSS 


Cx PROGRESSES AND SOLVES FOR THE PRESSURE FOR THAT MOLE 
C* 
180 Cem ye re COP Ra ae INC, RATIO, CRIT, CORR, 13 
NC = NC + 1 
IF (CABS (VAL) .GT.CABS(K( INC*N+1)#*2))THEN 
Crit = 5°CPir 
CORRN = CORN 
ip (eR DY Le. £ OSGeo To. 66 
Ve ie eo | 
GOL 4180 
ENDIF 
IF(CABS(CORR).GT.ERRHAX) GOTO 180 


1 Os@ re Gc) 
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PGs) (1G ar 


136 


20 
S17 @ 
of 


IF (CRIT. UP eC 
CRIT = C2 
GOTO 180 

ENDIE 

GOTO 57 

WRITE(4,S570)I 


IF (VAL.EQ.0) THEN 
ELSE 
CALL MODES2(VAL,DELZ,N,K, 2, CORR, INT, INC, RATIO, CRIT, CO-eae 
NGo= Ne ad 
IF(CABS(VAL) .GT.CABS(K( INC#'+ 1) B2) ) THEN 
CRIP =-5°CRI? 
CORRN = CORN 
Ie (CRIT. LT... 005) Co gemse 
VAL = Bloz()) 
GOTO 190 
ENDIE 
1F(CABS(CORR). GP ERaMax) come uceo 
IF (CRIT LT. CRIg2 a s. 
CRIT = CRITA 


GOTO. 190 
ENDIE 
GOTO 57 
WRITE(4,570)I1 
FORMAT(' NO CONVERGENCE! !!!! FIND THE ERROR’ .16} 
END IF 
EIG(1) = CSeRT( VAL) 
IDT? = IFIXM(DEPTHT/ (DEL 2 
Lome = 1 ee 
POIPE = P(IDT2)-~ PClDEL 
DDIEFS = DEPTET = (IDT1)~onE2 
PTK = P(IDT1) + PDIEF*DDISs, (gona 
IDR1 = IFIX(DEPTHR/(DELZ*2))*2 
IDR2 = IDR1 + 2 


PDIPF = PiIDR2 jue (iba 

= DEPTHR - (IDR1)*DEL 
PRX = P(IDR1) + PDIFE*DDITE. 
AIG(1) = CMELA(RESntete nie 
DENOM = INT*CSORT(AIG(1)*F1/2) 


CER KKK HRKRAKEKRKEKEKERERKRKRERR RARE RT RR TORE EES KT KKK RES eS RE Re ae 3 


C FOLLOWING INSTRUCTIONS OfbY FOR Cseis Ce ra Le lea Colt fies 

C DONE FOR RANGE OF CRITICAL ANGLES =D ECR DIFFERENT I@fPEDAl Ca 5 = ae 

CARRERE KRKEEK REET AT Re Ke eS ee ee eR ee Fe We ee 
ae CONTINVE 

C 

c 

Creare een kk eer kare xe eke eS RON RS OR oe ee, Re ee Ok, ee ee 

Ee THE FOLLOWING LOOP CALCULATES Thee Gtr 2A teil Cf Ta2esicy- eee 

GC DIFFERENT RANGES BY SUMMING THs COtie eee PRESSURE DUE ro Tre 

Ce INDIVIDUAL MODES 

Ce CAN BE SET UP TO PRODUCE ATTEiUe7 tO, eee 2CULAR RANco [iZERY 2 

Ce OR AT IRREGULAR INTERVALS “SUCE Ss SO et ea Ole Oe 

CF METERS, THEN EVERY 100 METERS TO i006 AND SO ON. 
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a~ 


c 
Cxxee 
Cx 

er 

C* 

Cc* 

ex 

c* 

Cr 

CcC* 


Ty ees ee ee 


IMAGECNELY(O. 1) 
DO 560 IRNG =1,1 


DO 561 NRNG = 300,300C0 
R=(10** IRNG) ~NRNG 
R = 50000 
PRESSR = O 
STRONG = 6 
DO 562 I = 1,N 
PRESS(!I)=PPRESS(1)*CEXP(IMAG*EIG(i)78)/SQRT(R) 
MAG(I)= CABS(PRESS(1!)) 
PRESSR= PRESSR + PRESS(1) 
PRG@y=GarG( PRESSS ) 
Pip(1) = AIMAG(FRESSR) 
FASE = ATAN(AIMAG(PRESSR) /REAL(PRESSR)) 
STRONG=CABS (PRESSR) 
WRITE (4,3434)1, PRESS(1),P2RESS(I) 
FORMAT(I6, 4215.7) 
PRbePo4e7) MAG(! ),StROts 
FORMAT (2515.5) 
Beh = se soemty ING jt 2m= (e1C(1)/COS(CRITA) )**2) 
REF = (RATEO) = FACT)/(RATIO + FACT) 
FRACT = CABS(REF) 


Come ee oe a (INGA jeez -bve(t)~*2)/( 2° k( 1NC*N+1)=SOTTC:M) 
OLOS = (1-FRACT)*CABS(COMP) 
PHAS = ATAN(AIMAG(REF)/REAL(REF) ) 


MoETA = —22e05(21G(1) /K(INC™N+1)} 


CONTINUE 
CONTINUE 
ATIEM = =20<°,LCC20(STRONS) 


CONTINUE 
CONTINUE 


CONTINUE 
3 LOOP 
END 


ie ee OE Sey ee Pease he peer wy lint ANC, RATIO,CAIT,CORRN, fy 
Selec e.  POoemmme ( LOCC] je yaa, COnS INT, E1,/A1,B1,F2,A2,E82,F3,A4 
iBo,52,8,A, WAG, AP, DENO 


RAKE RK TAKA RRR KKK HERRERA AEE KER wee ee RTH R TEKH RS ees 


GavowsUDROULPINE USES As eGea My ORDER RUNGE KUTIA TECHNIQUE 


ORC REGU Ar oe hoc ol noe ioe to DERIVATIVE SiLARTING aI THE SUREACE 
PDs eel iG hey eae DIlStancCe 2 TIMES TRAE DEPTH INCREMENT UNTIL 


GHESEOTTOM TS REACiaD 

Seer oon loin Pact Ok THE -Lr+1 STEP 

Chee oe ee ee eee temper No ooUn, UERIVATIVE AND INTECORSL 
Pesto eh Miao m PG Sine eeORKRSEC TIC TO Fe APPLIED TO THE 

SU] coo seo ole GuUloo> Of etde MORMGAL tTlODE HORI ZCHTAL WAVE- 
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CA Seas Scone iC Ween eve 

cA icoul ois TRE CORRECTION IS APPLIED AND [THE RESULTS1S -2Ss-73 
ce oO TSE “Mie = oes. 

CF IN ADDITICN THE PRESSURS AT ALL DEPTS INPSAVeLS Soe eee 

A THIS IS US#D IN THE CAbGCLAT ION Or Alas ee Ac 


C***=*SHOULD SAVE ONLY THOSE PRESSURES REQUIRED NOT; ALL INC*N+] 


Ct rt ence xb wa ek RS wR a a Bk ew on We 6 oe KR we A ee Nee 


PI = 33 feleso265¢ 


iNT = 070 
B= 0.0 
P( 1) = 35 
son tele 
DO” 18 hea NC eee 
Cc 
De = WELZ (Lei) 
Pl. = =B * (REL) 2 = 
Al = A + DELZ*EIL 
Bl = 3824 2e02-4 
C 
F2 = -B1*(K(u+1)- eee 
fh2 = Ao pe h2 eae 
B2 = B + DELZ “Al 
Cc 
ES = @82*(n(Let) 32 =e 
AS =e 2 Ee. 3 
BS = 5+ 2~PEna AZ 
Se 
F4 = =BS"(K(L+2)*"*2 = ieee 
B= B + 2°DELZ*(A + 2°h i 22 
P(L+1) =B 
A = A + 2*DELZ*(F1 + 2°62 275s + Faye 
IND =S5NT + B* "23 Pee 2 2 
cc 
18 COME Lave 
MiaAG = ChPUCA(O.G, 2.0) 
FILL = RESbI ht iNe=se1 ) jee 
DENC = CSORT( (FILL = VAL)/EEuL) 


AP = CSORT(EOL-VAL/COS (CREE? = s)he Ore 
CORR = (b= (47Ar 7 
VAG = VALe = COAR/COnG.! 
Al RETURN 
END 


m™ Ae 


SUBROLTINE MODESZ (VéAL,DELZ, i,',P,CORR, fit, DiC RATIC, Calo Cee 
COMPLEA K( 10001), ?( 10001) (Ve Oe eee eee ble eee ee 
1F3,A43,B83,84,5,4, HAG, veto 


Ces naw anes dee eeu Cee a re 
(ez THIS SUBROUTINE USES A FOURTH CREE Re RU Ge Ut Aw fee ao 

cc TO CALCULATE PRESSURE END ITS DERIVATIVE STARTING Sl sae 8s ae 
CF AND STEPPING UP A DISTANCE OF 2 Pitipse ee Geri INChkere 2 Cee 

Ce THE “SURPACE 1S ReEacoey 

Ce CALC AT EaCH L IS IN FACT FORT TRE ee 

Ce &T THE TOP, THE CALCULATED, PRESGURE Fee ast bye OND rece 

ce ARE USED TO DETERMINE THE CORRECT IOND SG gee hie teeta 
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ites ook one eos 

oe Soon Pewee ol eGUlooo Of PTH NORMAL »0DE HORIZONTAL WAVE- 

e* NUMBER Peewee ore oere pele eNO Lat RESULT 25S ralSeD lCAlK 
CF TO THE WaArtin PROGRAt 

aay ee eee ee eee otter Fl eG DEYTH: INTERVALS (S FPRESED EsCt 

ce Tae oe Oe ie eee CALCULATION OF ATTENJATION 

Crt eee Se exe KEAN RRR KR NR RRA NAR ANTS NASR RAS RN Te RR KR 


Pl ='3..2421592654 


INT = OQ 
Be BOOR— ee WwW +) 
B=) i 


P(INC*N) = B 


Ieee CHEE (0.0,1.0) 
Ellie] REAL(K(LELCOR) )*42 
BEN@e= CSORT( (FILL - VAL)/FILL) 


A = CSORT(FILL-VAL/COS(CRIT)**2)*B* IMAG/RATIO/DENO 
A = 0. 

Pore ct tM) Ay= 0. 

bons f= J. 1NC*N-1, 2 


c 
DP = D&LZ*(LFLOOR-L-2) 
Fl = -B * ({K(LFLOOR-&+1)**2 - VAL) 
Die= A Se DEbare 1 
Bl = B + DELZ*A 

c 
F2=0-61)*(K(25LOCR=L)**2 = VAL) 
LO8= Pee peL eae 
peo= Bebe ZA 

G 
PS = -62"(K(LELOOR-L)**2 = VAL) 
(ee A ee DEA 53 
So 1B ee DELZ*A2 

(fe . 
F4 = -83*(K(LELOOR-L-1)**2 - VAL) 
Pee ee Dee as(A + 2401 + 27a? + AD) /6 
PYLFLOCOR-b-2) = B 
Peete 27 DU LE (2) + 2-12 # 2°23 + FA) 6 
INT = INT + B**2*DELE*2 

e 

28 CONTINUE 


Sonne — S-A/ LT 

VAL = VAL + CORR,CORAII 
Zo RE2CRN 

END 


SUBROUTINE EIGENS(BOTTON, KOUNT,M,KAE, EIG2,RMIN) 
DIMENSION E(10001),D(10001),AK2(10001),EIG2(50),E1G(50) 
REAL KAE(10001) 


Cc PRINT <656, BOTTOM, MOUNT, M, RMIN ' 
CR KR RRR RR RE RE RE RE RR TERE RK RRR EE ERE RR 
cx m 
ee oo iO? Soe Ow ie Tels SUSROUZTINE BECAUSE IT CALLS IMSL oy 
Cx See. oer Geomeiiae noe LIiCAL VO MODIFY THAT FOR COMPLEY 
ee fiver) Gee lotNyeLUboeense [AME REFORE CALCULATED AND USED IN a 
Ce ODEs) AND MODESZ WHERE COMPLEX EIGENVALUES ARE USED mat 


oo 


4 

~ 
o 
r- 


: r Nisa a o~ . 
. i‘ Gee NUNe oe . . 
-_— = wae Sone ae & wo NOG, <x, wenticuierey 4 


(<x ee 
wa 

ha » 
a THIS SUERCUTINE 18 TO? Sees tee ee See een. = 


Cx > 


CR TRE wR we ke Cee RR AA NR RN Oe Te we ee OO OO OS OK Re ee Re oe eee ; 


D2 = BOTTOM 7 (KOUNT=1) 
DZ2 = "B25DzZ 


D(I) 
E(I) 
100 CONTINU 


0 
AR) 


C 
D(KOUNT) = 
x1) seoee 
TSH 
IER 
CALDME 


lod 


RT1LS (0,5,8OUNT Cia sh oe 


GO 6000 J = 1 fi 
BIGZ(J) = =DCJjye22 
Cree Re RRR EERE TT ew OR OT ee re ew OR OR Se OF Oe ee 
S eALTHOUGH THE FORTIULA IS EXAC? Loe SUBROUTINE BUNS IN SisGae 
CG PRECISION SO ON®& MUST BS (CARE EUL 
RSM KK RNR REO Re wR NEKO WEN we KE RRS RRR RRR ERKKR ARE eke 
EIG2EX = AK2(20) - (3.2471S92eSss-7e70- 1 yee SOrtcia 
6000 (COPE ius 
C 
ae leh Sei ©. 20) 
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or SDI 


COMPUTER PROGRAM FFP 


3] 
rg 
I 
o 
es, 

~) 
a: 
t— 


ro 


Pees: 


CO Perea bue hy, V1, DVI, K( 25000) ,EN,EIGVAL, VAL, 1MAG, FUNC, 
TAREA SIMPLE, SOUN, Vea, -EW, VRTVWN, WRONK 
REAL Ke SO), KAY(30), RATE(30),C(3C0) 


€ 
NUM=1 
PI = 3,14159265 
CMAX=0 
CMIN=2000 

CG 
FREQ = 50 


RMIN = 40 


€ INC , THE NUMBER OF INCREMENTS PER REAL MODE, MUST BE EVEN) 
CER RK RKKRKKK REA RKKR KEK RRR RR RK KE RARE NKR KARA RRR ARK KARR RR RA TN RTO 
ING = 2200 
12 READ (3,*) D(NUM),C(iuiM) 
Cc 


fee o NUM), . br. ©) GOTO 14 
KAY(NUM) = 2 * PI * FREQ / C(NUM) 
IF (C(NUM).GT.CMAM) THEN 

CMAX=C (NUM) 

DMAX=D(NUM) 
END IF | 
IF (C(NUM).LT.CMIN) THEM 

KMAX=KAY (NUM) 

DKMAX=D (NUM) 

CMIN = C(NUM) 

END IF 
IF (NUM .EO. 1) GOTO 13 


Beer 100 = 9 ( Ar (NUL je fay (eUM-1)) / (D(HUM)=D(TIUH-1)) 


13. NUM = NUM + 1 
Geno. 12 

14 BOTTOM=D(NUM-1) 
N = IFIX(2*BOTTOM*FREQ/CMAX+.5)- 
INCTOT = INC*N 
Ieee iNecren cl. 16000), 1NCTOT’ = 18000 
DELZ = BOTTOM/INCTOT 
IF (N.GT.50) THEN 

WRITE(4,100) N 


100 EORMAT( TOO MANY MOLES FOR THIS PROGRAM! ! MGREAL): =2 re) 
GOTO 11 
ENDS IEF 
ce 
CER RRR KR RK RK KK RR RR RK KAR ERR KERR KERR ERE RRS 
c 


Ome vVAGUE “Cr PAbSOR?r DION DEPENDENG ON FREQUENCY CNLY-NO TEMP DEPENDENCE) 
. 


CREEK RE KKKKRKKKKRKKR ERE RR HERE ERE RE KHER RK RK KR RRR RRR KEK R KR RK RK KKK KR TEES 


FFREQ = FREQ/1000 
ALFA1 = .0000092/(.7 + FFREQ**2) 
ALFA2 = .0046/(6000 + FFREQ**2) 


oa 


7} 
4 
tae 
{*) 
rc 
4 
Rm 
ry 
C 
7) 
Hy 
7 
le 
—_ 


ALFA3 = .QO000000+6 
ALFA = (ALFAL + AULPAZ * ALGas)* Sekee 2 
SLEa = .O0GOL 
G ALFA = 0 
C*x%**x*  EXAGGERATE ALFA TO SEE if HUMNSER Of DTITBReI CNS teca2 
NUM = 1 


DO 16 [C= Pea 
DP = (1G- 1) Dela 
K(IG) = KAY(NUM) + RATE(NUM) *(DP-D(NUZ1) ) 
LE (DP.GE..B CHUM + 1.) )- NUM "Senin ae 
iS CONTENVE 


e 
IF (DEPTHT.GT.DEPTHR) THEN 
DUP = DEPTHR 
DLW = DEPTHT 
ELSE 
DLW = DEPTHR 
DUP = DEPTET 
ENDIF 
é 
1Ul =| URIS (DUP 4 OsEZ<2) je 
[U2 = 10a = 2 
The TRIe (DENY (DELZ=2) \*2 
Th oe 
G 
C 
CRRA RARER RRA RRR KEEN HERRERA KARE REAR HK RAR TREK KR RARER RTT 
Ee ALFA IS CALCULATED FOR ONE FREQUENCY AND IS PART OF VAL(COMPLEX) 
cede Eke R EH RREHREE REE pew teens Ee ee 
E 
AREA = 0.0 
IFFP = 2047 
DO GG) =1) nee 
€ 
CK KEKE EWR AN RRR ETT RRR KAREENA eRe RK EK AA RRR AR RH ee 
G 
C DURING THIS PHASE OF THE PROGRAM MODES2 PROVIDES VALUES FCR 
e PRESSURE AND DERIVATIVE AT BOTH THE TRANSMITTER AND THE RECEIVER 
C STARTING WITH THE INITIAL CONDITIONS AT THE SURFACE 
c (IE PRESSURE = 0 AND DERIV IS AN ARBITRARY VALUE) 
C MODES 2 PROVIDES VALUES FOR PRESSURE AND DERIVATIVE FOR THE 
C LOWER POINT ONLY STARTING WITH THE INITIAL CONDITIONS AT THE 
é BOTTOM (IE DERIVATIVE = O AND PRESSURE IS AN ARBITRARY VALUES 
C A VALU FN REPRESENTS THE CONTRIBUTION OF EACH 'INCREMENT' 
€ AND THE AREA UNDER THE FN CURVE REPRESENTS THE TOTAL ‘PRESSURE 
Lg eed ddde kan aeanedbhenda aden tekken 
C 


R = 500 
WVNMCH = 1.1*KMAX/IFEP 

WVNM = 1*WVNMCH 

EIGVAL = CMPLX(WVNM, ALFA) 

VAL = EIGVAL**2 

CALL MODES1(VAL, DELZ, INCTOT, K, U,DU, IU], 1U2, DUP, DLW, V1, DV, 


a2 


ios oe salcees Gm aN Satie 
i thee, 

See hese (yn oro, DiCTOrT, 6 ¥y,/DV,I1L1,1iu2,DLW, INC) 

WRON- = V¥*bDVi=-V1-DV 
878 FORMAT (10E8.1) 
CER RRR TKKK RRR ERE RRR REE KKK REM KKK REE KER EKA KR KE RT KRE KR ERE HD NER TERK ATE TKR Ye eK 
C CALCULATE PARTIAL FIELD INTEGRAL FOR PARTICULAR VALUE OF 
et. WAVENUMBER COMPARE TO THAT CALCED BY PROGRAM 


GRRKK ERA KR AANA RARE ARAN AX KKKAK RHR AKA AK Ke KR RK RRA RAK RR TT 
CRivMem=weESORTIK(20)""2 = VAL) 

YEW = SIN(VRTWN*DUP) 

VEE = COS(VRTWN* (BOTTOM-DLW) ) 

WRONK = VRTWN*COS(VRTWN* B5OTTOM) 

SOLN = YEW * VEE / WRONK 

SIMPLE = U*V/WRON 

EN = U*V*CSORT(EIGVAL) /WRON 


OPC Or? (ie 


[MAG ="CMPLX(O,1) 

FUNC = FN*CEXP( IMAG*EIGVAL*2)*SORT(2/(PI*R)) 
CONT = CABS(FUNC 

AREA = AREA + FUNC*WVNMCH 


TOTAL = CABS(AREA) 
Ee iemi ion ae, Oe WRItE (2,190) °fI1GVAL,CONT, TOTAL 
TE (MOD (i eOe abe. 0) PRINT 190, EIGVAL, CONT, TOTAL 
55 CONTINUE 
EC 
11 STOP 
END 


RHA KKRTKKK TKK HR RK RRR NKR HK KA EKA RRR K KR RRR KR RETR REN KRKEUKKRKRK TKR ERRIKT 7 
VOeE Sieh OVIPeoe FESS AND eDeRIy 41 BOTH UPPER AND LOWER 
Pe UtS 


VAL - TRE INCREMENTAL HORIZONTAL WAVENUM - U2 TO 6192 VALUES 
Dee wee Bee ieee EN Oc (DEPIO, MUniB OF HODES, AND NUMNSEK 


OF INCREMENTS / MODE 

N - NUMBER OF MODES 

K - HORIZONTAL WAVENUMBER (COMPLEX) AT EACH DEPTH INCREMINT 
DUP - DEPTH OF UPPER OF TXNER AND RXER 

DLW - DEPTH OF LOWER OF TXER REER 

TU1/IU2 INCREMENT NUMBER FOR INC ABCVE AND BELOW 'UPPEK' 
ILI/IL2 INCREMENT NUMBER FOR It=!Cc ABOVE AND BELOW ‘LOWER’ 
INC - NUMBER OF INCREMENTS OF DEPTH 
OUTPUTS 

Ul - ‘PRESSURE’ AT UPPER 

DU1 - 'DERIVATIVE AT UPPER 

V1 - 'PRESSURE' AT LOWER 

DV1 - 'DERIVATIVE' AT LOWER 


CREA Ee eK EAE CRE Kate Ce Ek EC Ae Ae ea rAd RA Se Re ORK Rsv 7G we ew WO Oe Nay te oe ee 


COCNGY Go COO) Cat Gy Ca) 


SUBROUTINE MODES1( VAL, DELZ, INCTOT, K,U,DU, IU1, 1U2,DUP,DLVi,V1,DV1, 
Titel G2) 

COMPDEY K( 25000 9h (2 bs( 2) Bri(2),P1(2),F1,A1,B1,F2,A2,B2, 

Tao ese ae Ae PSIPe Plier GDIES, OIDIFF,U,V1,DU,DV1, VAL 


Ono 
-OOCO1 


B 
A 


oe 


18 


rap 


QAAANNA 


DO 18 {lb = were == 
Fl = -B « (N(L)#*2 - VAL) 
Al =: eee c= 2 
Bl = B + DELZ«*A 
PQ0= =Bi=(s(b-ljictee= oe 
A2 = A + DELZ * F2 
82 = B + DELZ *Al 
FS = =326(K(balj=°2)— oe 
A3 = A + 2*DELZ*53 
B3 = B + 2*DELZ*A2 
Fa = 353*(K(L+2)**2 Seen) 
5 = B + 2*DELZ*(A + 2*Al + 2=age Ali 7G 
A = A + 2*DELZ*(Fl + 2*F2 + 2°F3 + F4)/6 


TE((L+1).£Q.1U1) THEN 
poe 
peti) = s 
ENDIF 
IF ((L+1).EQ.1U2) THEN 
2 
Dei Zi. 6 
ENDIF 
IE ¢ (Geico tn)) rue 
PH) = 5 
Dei (ly aon 
ENDIF 
IF ((L+1).EQ.IL2) THEN 
Pie ae 
Douce 
ENDIF 
CONTINUE 
PDIFF = P(2) = P(1) 
PLIDIEE = Pi(2) — ely 
OODLES = DP(2) = eres 
OIDIEr = Dri]. oe 
DDIFT = DUP = (:U))}=Den- 
DIDIE: = DIS! = (uly pcr. 
U = P(l) + PDIFS*DDIFE/(DELa72) 
Vl = Pl(1) + PIDIFE*=DIDIFF/(DELZ*2Z) 
DU = DP(1) + CDIFE*DDIEE/( Danae) 
DV1 = DP1(1) + QIDIFE*DIDIFE/(DELZ*2) 
RETURN 
EN 


FOES eee ns ca, a Ra eg ees I MNase Aig et Late ee ee 


MODES2 PROVIDES PRESS AND DER AT 'LOWER' CNLY 
[Neuss 
VAL = THE IMCREMNENTAL HORIZONTAL WAVeN Ihe EO G1lo2 oc 


= 


DELZ- DEPTH INCREMENT EN OF DEPfH) Sis Cee Bes eee oe 


94 


ce) 
Oo 
G] 
< 
a 
tis 
755) 
© 
73 
2 
5 
}— 


& OF INCREMENTS, jeer 
Cc we NGNBER CS :.cDes 
€ eee eee ee Ube ae eemeGEX) AT S-Ci DEPTH INCREMENT 
C Pie br eieoe UePek OF THER AND RiER 
c Poe eo elem een Or Tin RAER 
eS ly eee eee USER FOR INC ASUVE AND BELOW ‘UPPER’ 
G Pee Pes CReie  MULISER FOR IMC ABOVE AND BELOW “LOWER ' 
c PNGe-UMEeR Ob 2NChEMENTS OF DEPTH 
o OUTPUTS 
Ss Ul - "PRESSURE! aT UP2ER 
Cc Buieee PERIVAT IVE 27 UPPER 
c V1 - 'PRESSURE' AT LOWER 
cS DV1 - 'DERIVATIVE' AT LOWER 
C 
Gwe te eee meee cK KKK KK KA KR EKER RK KEK RRA AKA ERR EEK KRARE RENT RRR KE KER EK 
C 
DeenoutINE MODESZ (Val, DELZ, INCTOT,&,V,DV,?ib1, tL2,DLW) 
Serie. K{(25000),2(2),DP(2),F1,A1,51,F2,A2,B2, 
ioe ess, 62,4), 8, PDIP ODIFE,V,DV,VAL 
B= .QO0001 
So 
BOe2e Ge — J INeror- 1, 2 
C 
EL = “ENG le fee =") 
Ba = Bee 2 = Ve) 
Al = Seu OE LZ 21 
Bl = 8 + DELZ=A 
c 
Pa >= -Bl*(R(Gle1)**%2 = VAL) 
a= Abt DELETE 2 
BZ. = B+ DELZ* Al 
C 
BS = )=B2*(k(bi-1)**2 = VAL 
Poe et eer 3 
Pome (> tee De AZ 
S 
Be OS = (Ri U1 =2)** 20> VAL) 
Ee So 29 pimeat(s = 2sAl, = 27A2 + 2354 ‘6 
eee eee ete 2 eS FE ) /6 
‘= 
Mala? yoree@. 151) TEN 
Bi =e 
DP(1l)= -A 
ae 
TE((L1-2).50.1L2) THEN 
EB 
DP(2) = -A 
EPDiS 
Ge 
28 Cor LIriue 
PDIFF = P(2) - P(1) 
ODIFF = DP(2) - DP(1) 
DORE ES=sDEW = (1L1)*DELZ 
ee ee OEE DoOlrE/ ( DELZ* 2) 
DV = DP(1) + QDIFE*DDIFF/(DELZ*2) 
Se 
2 RETURN 


END 


5) 


Or 
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